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Abstract: An occulter is a large diffracting screen which may be flown 
in conjunction with a telescope to image extrasolar planets. The edge is 
shaped to minimize the diffracted light in a region beyond the occulter, and 
a telescope may be placed in this dark shadow to view an extrasolar system 
with the starlight removed. Errors in position, orientation, and shape of 
the occulter will diffract additional light into this region, and a challenge 
of modeling an occulter system is to accurately and quickly model these 
effects. We present a fast method for the calculation of electric fields 
following an occulter, based on the concept of the boundary diffraction 
wave: the 2D structure of the occulter is reduced to a ID edge integral which 
directly incorporates the occulter shape, and which can be easily adjusted 
to include changes in occulter position and shape, as well as the effects of 
sources — such as exoplanets — which arrive off-axis to the occulter. The 
structure of a typical implementation of the algorithm is included. 

OCIS codes: (050.1940) Diffraction; (050.1970) Diffractive optics; (070.7345) Wave propaga- 
tion; (120.6085) Space instrumentation; (350.6090) Space optics. 
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1. Introduction 

One of the major goals of the coming decades is to directly image a terrestrial exoplanet. Di- 
rect imaging allows spectral information about the planet's atmosphere to be extracted from the 
spectrum of reflected light, including the potential presence of biomarkers, gases whose pres- 
ence and abundance may indicate the presence of life. Detecting planets directly is a difficult 
task, however, for two reasons: intensity ratio {contrast) and angular resolution. An Earth-twin 
which orbits a Sun-twin at 1AU emits 10 10 times less flux than the parent star in the visible 
band; if the system was located at 10 parsecs from Earth, the angular separation of the two 
objects would be lOOmas, which for most proposed space telescopes puts the separation under 
4A /D, Specialized methods are required to remove this flux at these small working angles. 

One proposed method of suppressing the starlight is an occulter, a spacecraft with a shaped 
edge flown in front of the telescope to block the starlight before it arrives at the telescope. 
The occulter size (tens of meters) and distance (tens of thousands of kilometers) are chosen 
so the angular extent of the occulter is smaller than some desired angle, on the order of 100 
milliarcseconds, so exoplanets in orbit about the star will still be visible. The edge is shaped 
to control diffraction, with the form chosen to suppress the light to a factor of 10 10 across the 
telescope aperture and over a wide spectral passband. 

A major advantage of this approach, compared to an internal coronagraph, is that it sig- 
nificantly loosens the tolerances on the optical quality and thermal stability of the telescope 
optics, as well as eliminating the need for a wavefront control system for active correction. 
Conversely, the occulter itself must be built, deployed, and flown in formation to its own set of 
tolerances, and demonstrating this capability is an area of active research Qj. These tolerances 
include limits on permanent errors in the occulter shape due to manufacturing and deployment 
errors; transient shape deformations from thermal fluctuations and vibration; and position and 
orientation changes of the occulter relative to the telescope-target axis from formation flying 

EH. 

The scale of an occulter — and the distance it must maintain with respect to its target — 
precludes optical testing on the ground. All of these effects must be modeled in order to build 
an error budget and verify that an occulter can satisfy these requirements. The capture of the 
large dynamic range in the resulting field — 10 s or greater in amplitude — drives the accuracy of 
the propagator used for the models, while the modeling of time-varying thermal shape defor- 
mations and closed-loop formation flying in broadband light drive the speed of the calculation. 
Some of the most in-depth simulated observations can take from hours to days to run, and the 
ability to quickly and accurately model propagation past an occulter is thus a major enabling 
factor in the characterization of occulter system performance. Section [2] gives an overview of 
the current techniques; the new technique presented in this work is described in Sections [3] 



and 3.1 and a suggested implementation scheme, along with a pseudocode overview of the 



algorithm and some computational results, are given in Section |4] 



2. Current occulter modeling techniques 



Occulter designs take the general form of a solid central disk surrounded by N p identical taper- 
ing structures {petals) around the edge, giving the whole structure the general appearance of 
a flower. (An example is shown in Fig. |T|) This shape is a result of their provenance as (0,1)- 
valued approximations to apodizers; the general method to design an occulter is to determine 
a smooth apodization profile A(r) which can provide the necessary starlight suppression, and 
then convert it into a binary shape with a sufficient number of petals so that performance is not 
undermined. The actual form of A(r) can be determined by optimization [4, 5| or selecting an 
existing functional form J6]|3. 



Fig. 1. An example of a typical occulter with 20 petals, with the associated Cartesian and 
polar coordinate systems. 

For an occulter designed in this manner, the region Q. in which the occulter is opaque can be 
written directly: 

Cl = {(r,e):0<r<R,0 G @(r)} (1) 
where 0(r) = 

n=0 

with (r, 6) being polar coordinates in the plane of the occulter. Each of the [...,...] is a set of 
points which defines the outline of a single petal; U denotes the union of these sets, which sets 
the boundary for the entire occulter. (We denote this boundary as <9£2, as it will be used later.) 



A(r), 1 A(r) 



(2) 



For a typical occulter, the occulter-telescope distance z = 5 x 10 7 m; the wavelength A = 
5 x 10~ 7 m; the occulter radius R = 25m; and the maximum excursions in the plane of the 
telescope and |tj| < 3m. These values place the occulter well within the Fresnel regime; 
that is, the region in which it is valid to approximate the exponent in the propagation integral 
with the first two terms of a power series (the Fresnel approximation). This approximation is 
suitable when JSJ. 

^»^K*-S) 2 +Cy-iJ) 2 ]L«- (3) 

Further, a standard propagation integral computes the electric field past a finite aperture; the 
occulter, on the other hand, is a finite obstruction with an effectively infinite aperture. Here, 
we can turn to Babinet's theorem [9|, which notes that field propagated through an aperture 
and the field propagated past its complement sum to the field that would have resulted were no 
obstruction present at all. 

If a plane wave with amplitude A and wavelength A is normally-incident on the occulter, 
as would be the case for an occulter properly aligned with a target star, the electric field at 
a distance z past the occulter would be written with the Fresnel approximation and Babinet's 
principle as 

U(4,17)=Aexp fl-^y^exp||i[(rcose-4) J + (rsin9-17) 2 ]Jrdrde') 

where £, and 77 are Cartesian coordinates at the downstream location. The "1— " term before the 
integral results from the use of Babinet's principle. Errors in (for example) occulter shape will 
modify i2, but Eq. (|4| will hold regardless. Computing the integral for points in the plane of a 
telescope aperture is the key step in modeling the performance of an occulter. 

One could attempt this integral by brute force: discretizing the shape of the occulter over 
a sufficiently fine grid to capture all of the necessary structure, and propagating this to the 
telescope aperture with direct integration or a matrix Fourier transform lflOl . This approach is 
not an easy task, as the occulter will be tens of meters across, with tolerances at the level of tens 
of microns[2], requiring a grid up to 10 6 x 10 6 or more to capture the full shape of the profile, 
not including array padding. 

A better approach is to take advantage of the structure of the occulter: in the absence of 
errors, it has an A^,-fold symmetry in the azimuthal direction, and most errors introduce small 
perturbations to that. Taking advantage of the symmetry leads to three main approaches to 
modeling propagation past an external occulter: 

1. Integrals in r: The Bessel function expansion [4| does the integral over explicitly, re- 
ducing the 2D integral to a series of single integrals in r whose contributions fall off ex- 
ponentially fast near the optical axis, which makes computation extremely quick in this 
region. The series can be modified to incorporate some shape errors [ 1 1 1 at the expense 
of slower convergence. Unfortunately, it is slower at locations far from the optical axis, 
and requires a good deal of effort to incorporate many types of errors, such as occulter 
tilt, into the model. Aside from the propagation modeling, the first, radially-independent 
term in this series is used in optimization approaches |0][5) to design occulter shapes. 

2. Integrals in : The Dubra-Ferrari integral [ 12 1 starts one step back from the Fresnel form 
of Eq. Q, at the Raleigh-Sommerfeld diffraction integral, and evaluates the integral over 
r to produce a pair of single integrals in angular coordinates, one of which holds for 
points inside the geometric extent of the starshade, and one which holds outside. This 



approach can incorporate most errors straightforwardly and spends the same amount of 
time to determine the field any point in the telescope aperture plane. Because of the 
manner in which the integrals are segmented, however, the integrands are more difficult 
to determine correctly for points physically outside the extent of the starshade. 

3. Perturbations to a nominal shape: The slit approach (131 assumes the electric field for an 
occulter with no errors has been calculated by one of the above methods, and approxi- 
mates the difference between the shape with and without errors as a series of long, thin 
boxes whose Fresnel transforms can be calculated quickly and exactly. This method is 
very fast but necessarily approximate. 

In this paper, we propose a fourth method — the boundary -diffraction-wave approach — which 
seems to out-perform the three extant methods for general purpose use. It maintains the 
constant-time property of the Dubra-Ferrari integral with a 2-3 x improvement in computa- 



tional speed; see Sec. 4.1 for further discussion. It incorporates all errors on the shape of the 
occulter natively, and can be written straightforwardly to include position and orientation er- 
rors, as well as model waves from off-axis sources, such as an exoplanetary system or the finite 
diameter of the star. It can also take advantage of the structure of the integral to evaluate the 
field at multiple wavelengths with little additional computation. The approach of defining the 
occulter solely in terms of its boundary is a good fit for tolerancing the sorts of errors that are 
expected for an occulter. Manufacturing errors can be applied by specifically deforming de- 
sired edge regions; errors in petal placement can be included by translating and rotating only 
the boundary points associated with that petal. 

3. The boundary diffraction wave formulation 

Most diffraction approaches to propagation through a finite aperture find their basis in the 
Huygens-Fresnel principle, which expresses the electric field at a point following the aperture 
as the sum over spherical waves emitting from every point within the aperture. However, a dif- 
ferent representation was introduced by Maggi and Rubinowicz which splits the field into two 
parts, a part which propagates according to geometrical optics and a part which incorporates 
the effect of diffraction from the boundary. This boundary-diffraction-wave (BDW) representa- 
tion was codified by the work of Miyamoto and Wolf 031 Q3), who showed that the boundary 
integral could be represented for any incident field entirely by a line integral around a vector 
potential. We will review their analysis briefly here. 

Miyamoto and Wolf lTT4l begin by writing the Kirchhoff formulation (see e.g. [9|) for the 
diffraction integral from the aperture. Consider a volume of space, bounded by a surface S; the 
field at any point P within the volume may be expressed as an integral over all surface points Q 
on S: 

Here we let s be the distance between P and Q, n be a vector normal to S at point Q, and Vg is 
the gradient evaluated at point Q. In the case of the occulter, the volume is the space z > 0, and 
S stretches over the z = plane, as well as having a component at infinity which is assumed to 
vanish. (Showing that this term vanishes requires some minor additional assumptions; see Sec. 
8.3.2 of (3). 

It can then be shown [ 14] that for V of this form, there exists a vector potential W, such that 



V = V x W, and for plane waves incident on an aperture, that associated vector potential is 



An s 1 +s ■ p 

U(P)= f[vxW(Q,P)-ndS. (7) 



s 

Here r is the vector from the origin to P; Q is a point within the aperture on S, and r' is the vector 
from the origin to Q. k = 2k /X is the wavenumber and A is the wavelength under consideration. 
The various s-variables become 

s = r'-r, (8) 

*=IM|, (9) 

s=*. (10) 
s 

and we assume the form of the plane wave in the vector direction p to be 

U (P)=Aexp(ikp-r), (11) 

with amplitude A. A coordinate diagram is shown in Fig. [2] 

Next, we consider for the moment the case of an opaque screen with an opening over the 
region £1, a subregion of S with boundary dQ.. We recall here that S is the entire bounding 
surface of the half-plane z > 0, including the z = plane; £2 is the section of it on which the 
occulter lies. (We are not considering the occulter case at the moment, but its complement; we 
will return to the occulter shortly.) In this case, the electric field from Eq. |7]) across the aperture 
is: 

U ap (P) = JJ VxW{Q,P)-nd£l, (12) 

Our approach is to reduce this integral to a line integral, using Stokes' Theorem. To do this, we 
consider singularities of the vector potential. 

In general, W will have singularities somewhere on S for a given P, as otherwise the field 
in the half-space past the aperture will be zero. In the specific case of the plane wave, the 
lone singularity occurs at 1 +s ■ p = 0, which occurs only when s is parallel to the propagation 
direction of the plane wave; thus, the singularity will fall in the aperture only for points P which 
fall in the beam as dictated by geometric optics. (Fig. [3] shows an example of this.) Miyamoto 
and Wolf |[T5l show that when Stokes' Theorem is applied, we get two distinct cases for the 
field at point P, depending on the singularity location: 

Uap(P) = A exp (ikp- r) + U^- B \P) for points inside the extent of the aperture (13) 
= U {B) (P) otherwise, 

where lj( B \P) is the counterclockwise line integral about the edge of £2: 

U (B) {P)= ( W(Q,P)-ldi. (14) 

J dQ. 

Here, £ is a unit vector in the direction tangent to the edge at any point, and d£ is a differential 
element of the boundary. 

This formulation gives the field following an aperture, but in our case we have the opposite: 
the occulter is opaque, and the surroundings are open. Here we return to Babinet's theorem. As 



we started with the plane wave in Eq. (Hi, this implies our occulter field U(P) satisfies 

U (P)=U(P)+U ap (P), (15) 



Fig. 2. A diagram of the points, coordinate systems, and vectors used in this paper. The left 
grid is in the plane of the occulter, and the right grid in the plane of the telescope aperture. 



or 



U (P) =— W(Q,P) • £d£ for points in the geometric shadow of the occulter, (16) 



This equation is the general form of the occulter field, in the boundary-diffraction-wave formu- 
lation. 

3.1. Vector potentials for occulter propagation 

From this suitably general representation, we can make appropriate substitutions for the specific 
case of an occulter; these will simplify computation of electric fields later. We define the vector 




The part of Eq. ( 16 1 which depends on geometric optics is termed (P): 



jj( g ) — o for points inside the extent of the aperture 
= Uo(P) otherwise. 



(17) 



Fig. 3. A diagram showing the projected extent of the aperture given by geometric optics. 
A plane wave incident on the aperture in plane A with vector direction p produces a shifted 
aperture in plane B. Two angles which parameterize the direction of the plane wave, 
and y/2, are shown as well. 



from the origin to P as r = , rj,z); r', the vector from the origin to Q, is defined as r' = (x,y,0). 
(See again Fig. [2]) We can write the s-variables explicitly as 

s = r'-r=(x-^y-r},-z), (18) 
*=||s|| = [z 2 + (x-§) 2 + (y-Tl) 2 ] 1 /2 (19) 

s= S . (20) 
s 

We consider first the case where the plane wave is normally-incident on the occulter, as 
starlight would be when the system is correctly aligned. Given this, our vector terms can be 
written as 

P= (0,0,1), (21) 
pr' = 0, (22) 

fxp=i(y-Ti,-(jc-S),0), (23) 



and the potential in Eq. (|6]i becomes 

exp (iks) = exp (ikz) exp [ik(s — z)] , 

W(Q,P) = —A exp (ikz f Xp[ik{S - Z)] S x p. 
\% s — z 



(25) 
(26) 



More generally, an off-axis source at an angle of y/\ off-axis and y/2 from the x-axis will have 

p= (siny/i cos y/2,siny/i siny/2,cosy/i), (27) 
p • r' = x sin y/\ cos y/2+y sin y/i sin y/2 . (28) 

Before the expression for the vector potential becomes too complicated, we will define three 
intermediate variables (/, g, and h) which we can substitute into Eq. (19i; these choices will 
greatly simplify subsequent computations. 



/ = zsiny/i + (x— cos y/2 cos Yi + (y — 1?) sin y/2 cos , 
g = -(x-§)sirn/A2 + 0-Tj)cosi^2, 

h = zcos y/i — (x — cos y/2 sin y/i — (y — 17) sin y/2 sin y/i = s p. 
Substituting these variables gives: 

s=[f 2 +g 2 + h 2 ] l/2 , 
h 

5 

s 

j'xp = - (/sin y/2 +g cos W2 cos Vi> 

— /cos 1/2 + g sin y/2 cos y/i , 
-g sin y/i) 

exp(ikh) — exp [z'fe cos y/\ — ik(x — <^)cos y/2siny/i — ifc(y — rj) sin y/2 sin i/^i ] 

= exp (— ikp ■ r') exp (ikzcos y/i)exp \jk{£, cosy/2 + f] sin y/2) siny/i], 

. . 1 , ., . . . . . exp \ik(s — h)] 

W(g,P) = — A exp (ifecos w ) exp \ik(c cos y/2 + T7 sin y/2) sin y/i - — - 

4/T .... v /; 



l+s-p= 1 
1 



(29) 
(30) 
(31) 

(32) 
(33) 

(34) 



(35) 



srxp, 



(36) 



a similar relation to Eq. (26 1. This equation is still exact; for y/\ <C 1, the usual case for objects 
in an exoplanetary system, we can note that 



1 _ 1 

s(s-h)~ p+g* 

f+g 2 



1 



[1 + ^]V», 



/ 2 +s 2! 



2// 



s — h 

and so the potential becomes 

W(2,P) s» — Aexp(fccosy/i)expte(i* cos y/2 + 77 sin y/2) sin y/i 

27T 

v(/,g) =(/ sin y/ 2 +^ cos y/ 2 cosy/i, 
— /cos y/2 + g sin y/2 cos y/j , 
-gsiny/j) = sxp. 



exp (iktyjp) 



f 2 +i 



(37) 
(38) 



v(/,g), 
(39) 
(40) 



This form proves to be particularly useful, as (f 2 + g 2 ) and v • i can be calculated without 
loss of precision, as neither will be of order z for the small y/\ case. A typical occulter might 
have z = 5 x 10 7 m, R = 25m, and || | and \ r\ \ < 3m, and k « 10 7 m _1 ; a representative exoplanet 
target for this occulter might have y/\ = 5 x 10~ 7 rad. We note that exponential terms of order kz 
have been separated from smaller terms; these should be calculated independently, as kz ~ 10 14 , 
and combining terms will lose precision in evaluating the exponent. 

4. Efficient implementation 

Evaluating the field at each downstream point requires the determination of two parts: the ge- 
ometric field and the boundary field. We assume, to begin, that we are given a set of N points 
representing the edge of the occulter which form a simple closed curve, as well as a list of M 
downstream points and a list of L wavelengths at which the field is to be determined. We also 
assume that the edge is sufficiently well-sampled that the section of the edge between explicitly- 
specified points is well-modeled by a linear segment, as this allows us to use the midpoint rule 
for numerically approximating the integral. Given the slowly-changing shape over the majority 
of the petal, this is a good assumption; small regions such as petal tips may be specified more 
finely. 



The geometric field is still a plane wave of form Eq. (Hi outside the aperture, and is 
everywhere inside; finding the geometric field thus becomes a problem of checking whether 
a downstream point falls behind it or not. We do this by treating the edge as a polygon and 
finding the winding number of the downstream point (<!jo, Tjo); this can be done efficiently with 
an 0(N) routine such as polywind in Numerical Recipes lfl6l . and holds regardless of the 
complexity of the occulter shape. This approach also speeds up multiband calculations, as the 
geometric extent of the shadow determined this way is wavelength-independent. Lateral errors 
in the occulter-telescope position may be included by adding a At; and Arj to all points at 
the telescope plane, and that should be done first. Note that in the case of an off-axis source, 
the geometric shadow is shifted laterally as well; in this case, the winding number should be 
calculated around (<i;o — zsin cos y/2 7 ??o — zsin y/j sin y/2). We then have: 

U ( G) (<jjo > i?o , z ) = A exp (ikz cos yi ) exp [ik ( £ cos y/ 2 + rj sin y/ 2 ) sin y/i ] , (41) 
for points outside the shadow 
= otherwise. 

In some cases — for example, when the telescope aperture remains close to the center of the 
occulter and the perturbations to the ideal shape of the occulter are small — the winding number 
may be able to be determined independently. This determination should be done if possible, as 
even an efficient algorithm can increase runtime significantly for boundaries containing a large 
number of points. 

The boundary field can be evaluated efficiently with midpoint-rule quadrature running around 
the edge. First, for each of the N line segments running counterclockwise around the occulter 
edge, we derive the midpoint of line segment j going from (xj,yj,Q) to (xj+i,yj+\,0): 

{xf\yf\o) = ([xj+x j+l }/2,\yj+y j+1 ]/2,0), (42) 
and the vector pointing along the segment 

(xf,yf,0)j = -Xj,y j+1 -yj,0). (43) 



with (xj + i,yj + i) equal to (xi,yi) when j = N. This derivation need only be done once at the 
beginning, as it holds for all M points and L wavelengths. (We note that, while the midpoint rule 



has been used for the quadrature, this is certainly not a requirement, and higher order methods 
could be used; Eq. ( |39| ) will still hold. Using a higher-order method may come at the cost of 
runtime.) 

For each point in the downstream field, we calculate /, g, and h for all N segments following 
Eq. pj: 



fj = z sin + (x) - 4o) cos y/ 2 cos y/i + (yj 
gj = - {x [ f - sin y/ 2 + - i?o) cos y/ 2 
/;,■ = zcos v/i - (.^ m) - £0) cos y/2 sin y/i - (y) 
and we then need only two intermediate terms: 

T tj = fj+S% 



(m) 



Tjo J sin v/2 cos y/i, 



770 ) sin y/2 sin y/i , 



(44) 
(45) 
(46) 



(47) 



T 2j = y(fj,gj).iM 



y (z sin y/i sin y/ 2 + (jf' ~ ^o) cos y/i ) - yf (z sin y/i cos y/ 2 + (x^ - £ ) cos y/i ) , 



(48) 



to give us: 



U (B) (^0 , *7o , z) = —A exp (ikz cos w ) exp [/fc(| cos y/ 2 + r/ sin y/ 2 ) sin m] 
2K 



n exp (ifc^) 

x L 

7=1 



71; 



(49) 



Note that all of the steps prior to Eq. ( 49 1 are wavelength-independent; different values of k 
may be iterated over in the final step. It is advisable to keep terms of order kz (e.g. fecos y/\ 
or kh) separate when evaluating to maintain numerical precision. A typical implementation 
might take the form shown in Algorithm[T] The for-loops over j in particular are well-suited to 
vectorization where this is supported. 



for j = I to N (all occulter edge poin ts) d o 

I Find Xj , yj , xj' , yj' using Eq. ( 42 1 and Eq. ( 43 1 
end 

for q = 1 to M (all points at telescope aperture) do 
Apply lateral error to point q, if any 
Get winding number around point q with polywind 
for j = 1 to N (all occulter edge poi nts) do 



I Find fj,gj,h j ,T 1 j,T 2 j using Eq. (44 1 through Eq. (48 1 
end 

for s = 1 toL (all wavelengths) do 

Find f/' B ) by summing over hj,T\j, T 2 j in Eq. \49 \ 
Find U {G) using the winding number and Eq. (Wlb 
Find U = — for point q and wavelength s 
end 



end 



Algorithm 1: Pseudocode representation of BDW calculation 



4.1. Computational results 

In the standard formulation of the Dubra-Ferrari integral (DF), points outside the geometric 
shadow are calculated with a different set of integrals which require searches for nearest- and 
further-neighbor points along the boundary for each point at the telescope aperture. To simplify 
direct comparison between the Dubra-Ferrari approach and the boundary-diffraction-wave ap- 
proach for this paper, we choose to select a region which lies entirely in the geometric shadow 
of the occulter. This selection also fixes the winding number to 1 at all points under consider- 
ation, letting us bypass an explicit calculation. This choice is not an unreasonable assumption; 
the telescope will generally be entirely in the geometric shadow for a properly aligned occulter- 
telescope system. 

Both approaches have similar formulations, with the primary difference coming in the fact 
that the Dubra-Ferrari integral uses angular coordinates 9 centered about (§0)f?o)> an d so tne 
angular distance between occulter edge points must be rederived for each point at the telescope 
aperture. This calculation turns out to be one of the major sources of overhead during the DF 
calculation. 

Propagation runtime for a 128x128 output field 
10 3 , — — — — — u 




Number of boundary points 



Fig. 4. A comparison of the time to run the two propagation algorithms, as a function of 
number of points along the edge of the entrance aperture. 

To compare the runtimes of the Dubra-Ferrari and BDW implementations, we consider a 
closed occulter boundary specified by 164000 individual points, corresponding to the occulter 
shown in Fig. [T] We select subsets of these points, and run Matlab implementations of both 
algorithms on this subset. Fig. [4] shows the time to implement the propagation as a function 
of number of boundary points; the BDW algorithm remains 2-3 x faster over a wide range 
of polygon sizes. Fig. [5] shows the intensity profile in a region around the telescope aperture, 
plotted on a base- 10 logarithmic scale. The plots of this intensity profile created with Dubra- 
Ferrari and BDW are not visually distinguishable; within the telescope aperture, the maximum 
difference in intensity is 4.8 x 10~ 15 . 
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Fig. 5. Intensity at the plane of the telescope aperture following the occulter, plotted on a 
base- 10 logarithmic scale. The telescope aperture is shown by a dashed circle. 

5. Conclusion 

We present a new approach to modeling fields following planet-finding occulters, based on 
the boundary-diffraction-wave formulation laid out by Miyamoto and Wolf. It has proved to 
be simple to modify the shape and recalculate with this approach to include various types of 
occulter errors, as well as fast to evaluate, improving on the current best method by 2-3 x with 
no loss in precision. An efficient implementation has been laid out, and we hope it proves useful 
to others in the field. 
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